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ABSTRACT 

This paper considers the phenomenon of deep core collapse in collisional stellar sys- 
tems, with stars of equal mass. The collapse takes place on some multiple, of 
the central relaxation time, and produces a density profile in which p oc r~ Q , where 
a is a constant. The parameters a and £ have usually been determined from simpli- 
fied models, such as gas and Fokker-Planck models, often with the simplification of 
isotropy. Here we determine the parameters directly from iV-body simulations carried 
out using the newly completed GRAPE-6. 
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1 INTRODUCTION 

Consider a spherical non-rotating stellar system in dynamic 
equilibrium. Two-body encounters drive a slow evolution 
of the system. In response to the well known gravothermal 
instability of such systems (Antonov 1962, Lynden-Bell & 
Wood 1968, Hachisu & Sugimoto 1978) the core contracts. 
Eventually the central relaxation time is so short that the 
■ core loses thermal contact with the outer parts of the sys- 
tem. Thereafter the central parts of the system evolve in 
a self-similar manner, unaffected by boundary conditions 
(Lynden-Bell & Eggleton 1980). 

In this self-similar regime all central parameters evolve 
as powers in r, where r is the time remaining until collapse 
ends. Here we neglect variations in the Coulomb logarithm 
in the expression for the relaxation time. If t rc is the central 
relaxation time (defined as in Spitzer 1987, eq. 2-62), and 

Pc 

Pc is the central density, it follows that — t rc — f, where £ is 

pc 

constant. At the same time the density profile approaches 
a power law p oc r~ a , where a is another constant. 

Very little can be said about a and £ on general grounds. 
Since the core is nearly isothermal we may expect that the 
evolution timescale is much larger than t rc , and so £ << 1. 
Lynden-Bell & Eggleton showed on physical grounds that 

2 < a < 2.5. It has been claimed (Lancellotti & Kiessling 
2001) that a = 3, on the basis of the scale invariance of the 
Fokker-Planck equation. It is shown in appendix A that this 
is too restrictive a condition. 

Precise determination of the parameters a and £ have 
been obtained by a variety of methods (Table 1). Lynden- 
Bell & Eggleton themselves determined parameters equiva- 
lent to a and £, using an isotropic gaseous model of a stel- 
lar system, by determining the self-similar solution directly. 



This is an eigenvalue problem in which their two parame- 
ters are eigenvalues. In common with all gaseous models, 
the result for £ depends on a constant which is usually de- 
termined by comparison with results of some other method, 
and so the value is not given in the Table. Much earlier, Lar- 
son (1970) determined equivalent parameters by analysing 
a time-dependent solution of an anisotropic model based on 
moments of the Fokker-Planck equation. 

These two methods (eigenvalue problems and analysis 
of time-dependent solutions) have been applied by a num- 
ber of authors using various models and are listed in Ta- 
ble 1. Where necessary their results have been converted 
to yield values of a and £ using relations among core pa- 
rameters developed by Lynden-Bell & Eggleton, except that 
we uniformly use Spitzer's relaxation time. For the theo- 
retical models, the value of the Coulomb logarithm used 
in the model is assumed to exactly cancel that in Spitzer's 
definition. In several cases no values were given by the au- 
thors themselves, and so we have added notes to indicate 
our source for the values given. Only systems with stars of 
equal mass are considered here. 

While earlier discussions of this topic (e.g. Spitzer 1987, 
Louis 1990, Louis & Spurzem 1991) were restricted to com- 
parison of results from simplified models, our main aim in 
this paper is to add data from new A-body simulations. 
These and earlier A-body results are identified in col. 4 of 
Table 1 by the value of N used, or the range of N. 



2 RESULTS OF A-BODY SIMULATIONS 

We have performed a new series of A-body simulations of 
isolated clusters starting from Plummer profiles and contain- 
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Table 1. Determinations of a and £ 



Source 


a 








Model 


Larson (1970) 


2 


41 1 





00160 2 


anisotropic moment 


Louis (1990) 3 


2 


20 





00212 


isotropic moment (eigenvalues) 


Louis (1990) 3 


2 


23 





00123 


anisotropic moment (eigenvalues) 


Lynden-Bell & Eggleton (1980) 


2 


208 






isotropic gas (eigenvalue) 


Louis & Spurzem (1991) 3 


2 


23 






anisotropic gas (eigenvalues) 


Cohn (1980) 


2 


23 





0036 


isotropic Fokkcr-Planck 


Heggie & Stevenson (1988) 


2 


23 





00364 


isotropic Fokker-Planck (eigenvalues) 


Takahashi (1993) 


2 


23 





00365 4 


isotropic Fokkcr-Planck (eigenvalues) 


Cohn (1979) 


2 


27 





006 


anisotropic Fokkcr-Planck 


Takahashi (1995) 


2 


23 





0029 


anisotropic Fokker-Planck 


Duncan & Shapiro (1982) 


2 


2 





006 5 


Monte Carlo anisotropic Fokkcr Planck 


Joshi ct al. (2000) 


2 


2 






Monte Carlo anisotropic Fokker Planck 


Gicrsz & Heggie (1994) 


2 


17 6 






N = 500 (average of ~ 50 cases) 


Makino (1996) 


2 


36 






N = 32k 7 


This paper 


2 


26 





0030 


N = 8k - 64k 



Notes 

1. From eq. (43), p c 

2. From cq.(44) 

3. Several variants are considered in these papers. The values quoted are those highlighted by the authors in the abstract 

4. Assumes relation for / (distribution function), p c and central velocity dispersion for a Maxwcllian. 

5. Depends on the assumed value of the Coulomb logarithm; sec Spitzer (1987, p. 95) 

6. From Fig. 14 

7. From Fig. 4 



ing N = 8192 to 65536 stars. The simulations were carried 
out on the recently finished GRAPE-6 boards at Tokyo Uni- 
versity, using a specially adapted version of the fully colli- 
sional A-body code NBODY4 (Aarseth 1999). All runs were 
performed well into the post-collapse phase. Details of the 
runs can be found in Table 2. 

We first determined the position of the cluster center, 
using the method of Casertano & Hut (1985). According to 
Spitzer (1987, p. 149), the average density p c inside the core 
radius r c is 0.517 times the central density p(0) in an isother- 
mal model. Since also the following relation, connecting r c , 
the 3d- velocity dispersion in the core a c , and the central 
density, holds for an isothermal model: 



°l = yGp(0)rc , 



(1) 



the core radius can be estimated by the following relation: 



0.517 



3< 

4 TlGpc 



(2) 



Although our clusters initially do not follow isothermal den- 
sity profiles, their cores approach isothermal profiles as core- 
collapse proceeds. To calculate r c , we first estimated the 
average density p c inside the core and the 3d-velocity dis- 
persion from the innermost 1% of the stars and obtained a 
value for r c from eq. 2. a c and p c were then calculated using 
the stars inside r c and the whole procedure was repeated for 
10 iterations. We found that this was enough to determine 
r c and a c with sufficient accuracy, since after 10 iterations 
the relative changes of these values between successive itera- 
tions were of order 1% or less. Finally, the core collapse rate 
£ was calculated by comparing the core density at two suffi- 
ciently separated points in time and calculating the central 
relaxation time at the midpoint. Fig. 1 shows the evolution 
of £ as a function of the scaled energy xo = -30(0) /al for 




Figure 1. Core collapse rate £ as a function of central escape 
energy xo for the 64K run. Points denote individual values cal- 
culated for all times when data was stored. The solid line shows 
the run of the mean log(£). £ decreases steadily as core collapse 
proceeds, and becomes constant at around xq = 10. 



the 64K model. Here 0(0) is the central potential, calculated 
by excluding the star nearest to the cluster centre. 

The A-body data is noisy since the time derivative of 
the central density is used for calculating £. Nevertheless, it 
can clearly be seen that £ decreases until about xo = 10, af- 
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Table 2. Results for the core-collapse time Tqc, a an d £ from 
JV-body simulations 



N 


8192 


16384 


32768 


65536 


N Sim 


8 


3 


2 


1 


< T CC > 


1967 


3640 


6796 


12218 


<a> 


2.24 


2.26 


2.28 


2.26 


<a> 


0.0030 


0.0031 


0.0030 


0.0029 



ter which it becomes nearly constant. This is in good agree- 
ment with the behavior found by Cohn (1980) and Takahashi 
(1995) in their Fokker-Planck simulations. Core-collapse is 
completed at around xo = 13.0. Taking the mean £ from 
all data points with xo > 11, we obtain a limiting value of 
£ = 0.0029 for this run. Similar values are obtained for runs 
with other N (see Table 2). Taking the mean over all per- 
formed runs, we obtain £ = 0.0030, which is in good agree- 
ment with what Takahashi (1995) obtained from anisotropic 
Fokker-Planck calculations (see Table 1). 

In order to measure the density gradient at the time of 
core- collapse, we proceeded in the following way: The time 
of maximum core contraction was determined from the time 
when the potential energy at the cluster centre (4>{0), cal- 
culated as above) reached its first minimum. We then cal- 
culated the stellar density as a function of distance from 
the centre and fitted power-law distributions to the den- 
sity profile inside 0.1 half-mass radii. The slope of the best 
fitting power-law was determined by a KS-test. Mean val- 
ues of a can be found in Table 2. In general, we obtain 
somewhat larger values than the Fokker-Planck calculations, 
and our results seem to be only marginally compatible with 
the a preferred by most Fokker-Planck and gas calculations: 
a = 2.23. For the range of particle numbers studied, no clear 
change of a with N can be seen. 

Fig. 2 compares the combined iV-body data from runs 
with N = 8192 to 65536 stars with various power-law pro- 
files. In order to determine the slope of the density profile, we 
fitted data up to a maximum radius of r — 0.1 rnaif- Fig. 2 
shows that outside this radius, the density profile cannot 
be fitted by a single power-law any more, while inside from 
there the value of a will not depend on the maximum radius 
used for the fit. Inside r = 0.1r Ha if, we obtain a slope of 
a = 2.26 for the density profile. A value of a = 2.23 would 
give a bad fit to the combined iV-body data, but might be 
possible in the N — > oo limit if a is changing slowly with N. 
A value of a = 3.0 is completely ruled out by our iV-body 
simulations. 

Simplified anisotropic models for star cluster evolution 
(e.g. Giersz & Spurzem 1994) predict that, at a given radius 
inside the self-similar regime, the anisotropy approaches a 
finite value at the end of core collapse. Our results for the 
evolution of the anisotropy profile are shown in Fig. 3. At 
the end of core collapse, we obtain an anisotropy profile 
closely resembling Fig. 2 in Giersz & Spurzem, though the 
value of the anisotropy parameter A at small radii is a little 
smaller, around 0.2. The maximum value of about 1.1 (which 
occurs outside the self-similar regime) is a little larger than 
theirs. The time evolution of the anisotropy within different 
Lagrangian shells closely resembles the results from averaged 
1000-body models shown in their Fig. 11. 



0.4 

p(r) = r " with a = 2.23 

a = 2.26 

a = 2.29 

0.2 




Figure 2. Density profile of the combined set of iV-body runs 
compared with 3 different power-laws p(r) = r~ a . All data is 
divided by the density of a power-law profile with p(r) = r~ 2 - 26 . 
Power-law profiles are adjusted such to contain the same number 
of stars as the iV-body data inside r = 0.1 r Half - The best fit 
is achieved for a = 2.26. The uncertainty of this value does not 
seem to be larger than 0.02, ruling out a value of a = 2.23 in the 
10 4 < N < 10 5 range. 







"1 

T 


'1 1 

= 







T 


= 11000 






T 


= 12200 






T 


= 12218 //" 








// 


\ \ 






// 


X \ 






// 


\ \ 






// 


\ 






// 


\ 






// 




\ _ 




// 
















J 1 








/ 

/ 








/ 








0.01 0.1 1 10 100 

r [NBODY] 



Figure 3. Anisotropy profile for the run starting with N = 65536 
stars for 4 different times. The anisotropy parameter A rises as 
core-collapse proceeds until it becomes nearly constant inside the 
self-similar regime. 
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3 CONCLUSIONS 

Late core collapse of stellar systems is one of the few regimes 
where the evolution becomes relatively simple. For systems 
with particles of equal mass, simple arguments imply that 
many properties (central density and velocity dispersion, 
core radius, density profile, etc.) approach simple limiting 
forms, which can be characterised by just two parameters. 
A common choice for these is the dimensionless rate of in- 
crease of the central density, £, and the index of the power- 
law dependence of density on radius outside the core, —a. 

In this paper we explain why a is not determined, as 
has been argued, by the scale invariance of the Fokker- 
Planck equation. We review historical determinations of 
these parameters based on numerical solutions of the Fokkcr- 
Planck equation and other simplified models for the evolu- 
tion of stellar systems. Our main contribution, however, is 
to present determinations of these parameters directly from 
new large direct iV-body computations. We find that the 
core collapse rate £ (= 0.0030) agrees satisfactorily (within 
the statistical error of the N-body results) with those deter- 
mined by better simplified models. The density profile index 
a, is slightly steeper, our best value being about 2.26. 
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APPENDIX A: SCALING OF THE 
FOKKER-PLANCK EQUATION 

It has been shown by Lancellotti & Kiessling (2001) that 
the Fokker-Planck equation admits a unique scale invari- 
ance, and that therefore the self-similar solution requires a 
limiting density profile p oc r -3 , i.e. a = 3. Here it is shown 
why the value of a is not determined by the scale invariance 
of the Fokker-Planck equation. 

Consider first the model problem 

2L+V^L f (A1) 

dt +V 'dr r3-dv _/ ' {A1) 

which can be interpreted as a Fokker-Planck equation for 
a distribution of Keplerian oscillators with energy E = 

-v 2 — 1/r. The right side of eq. (Al) is a simple collision 
term, chosen only so that it has the same scaling property 
as in the Fokker-Planck equation of collisional stellar dy- 
namics. Indeed eq. (Al) admits the unique scaling / — > fxf, 
t -> p~ 1 t, r -^i~ 2/3 r, v -^ 1/3 v. It follows that eq. (Al) 
admits self-similar solutions of the form 

/ = t~ 1 F (rt- 2/3 ,vt 1/3 ) , (A2) 

where F is some function satisfying a certain partial 
differential equation. Hence the space density is p = 
t~ 2 f F (rt~ 2/:i , v*) d 3 v*, where v* = vi 1/3 ' This is sta- 
tionary at large r only if p oc r~ 3 , i.e. a = 3. 

These are not the only self-similar solutions of eq. (Al), 
however. There are also solutions of the form 

/ = t~ l F (t f3 E) , (A3) 



where f3 is any constant and F satisfies a certain ordi- 
nary differential equation. Such solutions are certainly self- 
similar, in the sense that the function / evolves by time- 
dependent scalings of /, v and r. The reason why eq. (Al) 
admits a wider class of self-similar solutions than those of 
the form (A2) is that solutions of the form (A3) also satisfy 
the differential equation 

dr r 3 dv 

and hence also the simpler Fokker-Planck equation 

This pair of differential equations admits a much wider class 
of scalings 

f -> Vf, t -> p~ l t, v -s- vv, r -> v~ 2 r . 

Another interpretation of this situation is to observe 
that (A4) is the orbit-averaged version of eq. (Al) (cf. 
Spitzer 1987). In the language of stellar dynamics it is 
the equation obeyed by solutions which evolve slowly, on 
a timescale t ev much longer than the orbital or crossing 
timescale, t cr . Indeed the first term on the left of eq. (Al), 
and the term on the right, are of order f/t ev , while the re- 
maining terms of eq. (Al) are of order f/t cr . If we restrict 
attention to self-similar solutions of the form of eq. (A2) we 
are, in effect, insisting that t ev oc t cr . Indeed the scaling of 
eq.(Al) is exactly the same as that of the iV-body equations 




and the analogues of eq.(A2) are then the familiar homo- 
thetic solutions in which r, oc t 2 ^ :i (cf. Arnold et al. 1997, 
p.65). 

For the Fokker-Planck equation of stellar dynamics the 
situation is a little more complicated. Its general and orbit- 
averaged forms are 

and 

f+i<f >=<(!)> <-) 

where ( ^f-) is the collision term, d> is the potential, 
\otJ c 

E = <f> + -v 2 and () denote an orbit average. Eq. (A6) is 

not obtained from eq. (A5) by assuming that f = f (E, i) . 
Instead, we must assume that a solution of eq. (A5) may be 
expanded in the form / = /o + e/i + ... where e = t cr /t ev 
and fo = fo (E,t) . Then fo obeys eq. (A6), while /i and 
higher terms obey appropriate linearised forms of eq. (A5). 

For self-similar core collapse, fo may be taken to be of 
the form of eq. (A3). There is no reason to suppose that /i 
enjoys the same self-similar evolution as fo. Thus the solu- 
tion of eq. (A5) for core collapse is nearly self-similar (up to 
terms of order t cr /t ev ) but is not confined by the highly re- 
strictive scaling properties of eq. (A5) itself. As core collapse 
comes to an end, with a small number of stars remaining 
in the core, t ev decreases and becomes nearly comparable 
with t cr - Then new phenomena beyond the Fokker-Planck 
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equation become important, such as formation of binaries in 
3-body encounters. 
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